use rest_tensors::BasicMatrix;
use tensors::{MatrixFullSlice, MatrixFull};

use super::DMatrix;

pub enum C2S {
    L0(DMatrix1x1),
    L1(DMatrix3x3),
    L2(DMatrix6x5),
    L3(DMatrix10x7),
    L4(DMatrix15x9),
    L5(DMatrix21x11),
}

pub type DMatrix1x1 = DMatrix<1>;
pub type DMatrix3x3 = DMatrix<9>;
pub type DMatrix6x5 = DMatrix<30>;
pub type DMatrix10x7 = DMatrix<70>;
pub type DMatrix15x9 = DMatrix<135>;
pub type DMatrix21x11 = DMatrix<231>;
//pub struct DMatrix1x1 {
//    size: [usize;2],
//    indicing: [usize;2],
//    data: [f64;1]
//}

//Transformation from normalized Cartesian functions, to normalized spherical hamonic functions
// H. BERNHARD SCHLEGEL, MLCHAEL J. FRISCH, Transformation Between Cartesian and Pure Spherical Hamonic Gaussians
// International Journal of Quantum Chemistry, Volume 54, 83-87 (1995)
// From Table 1
// v(l,m)=1/sqrt(2)(r(l,m)+ir(l,-m)) => r(l,m) + i r(l,-m) = sqrt(2) v(l,m), therefore
// The real part of v(l,m) times sqrt(2) is r(l,m)
// The image part of v(l,m) times -i sqrt(2) is r(l,-m)
pub const C2S_L0: C2S = C2S::L0(DMatrix1x1{
    size: [1,1],
    indicing: [1,1],
    data: [1.0]
});

//pub struct DMatrix3x3 {
//    size: [usize;2],
//    indicing: [usize;2],
//    data: [f64;9]
//}

pub const C2S_L1: C2S = C2S::L1(DMatrix3x3{
    size: [3,3], 
    indicing: [1,3], 
    data: [
           //  x    y    z
               1.0, 0.0, 0.0,    // (1,-1)
               0.0, 1.0, 0.0,    // (1, 0)
               0.0, 0.0, 1.0     // (1, 1)
]});


//pub struct DMatrix6x5 {
//    size: [usize;2],
//    indicing: [usize;2],
//    data: [f64;30]
//}

/// The transofrmation matrix from Cartesian to spheric for L=2 (D)
///        let a1 =  0.866025403784438647; //sqrt(3)/2 
///        let a2 = -0.866025403784438647; //-sqrt(3)/2
///        let rel_max = [
///        //  0    1    2    3    4    5
///        //  x^2, xy,  xz,  y^2, yz,  z^2  
///        //  200 110  101  020  011  002
///            0.0, 1.0, 0.0, 0.0, 0.0, 0.0,  //(2,-2)
///            0.0, 0.0, 0.0, 0.0, 1.0, 0.0,  //(2,-1)
///           -0.5, 0.0, 0.0,-0.5, 0.0, 1.0,  //(2, 0)
///            0.0, 0.0, 1.0, 0.0, 0.0, 0.0,  //(2, 1)
///             a1, 0.0, 0.0,  a2, 0.0, 0.0,  //(2, 2)
///            ].to_vec();
pub const C2S_L2: C2S = C2S::L2(DMatrix6x5 {
    size: [6,5], 
    indicing: [1,6], 
    data: [
        //  0                     1    2                      3    4    5
        //  x^2,                  xy,  xz,                    y^2, yz,  z^2  
            0.0,                  1.0, 0.0,                   0.0, 0.0, 0.0,  //(2,-2)
            0.0,                  0.0, 0.0,                   0.0, 1.0, 0.0,  //(2,-1)
           -0.5,                  0.0, 0.0,                  -0.5, 0.0, 1.0,  //(2, 0)
            0.0,                  0.0, 1.0,                   0.0, 0.0, 0.0,  //(2, 1)
            0.866025403784438647, 0.0, 0.0,-0.8660254037844386472, 0.0, 0.0,  //(2, 2)
]});
//pub struct DMatrix10x7 {
//    size: [usize;2],
//    indicing: [usize;2],
//    data: [f64;70]
//}

/// The transofrmation matrix from Cartesian to spheric for L=3 (F)
/// let b1  =  1.060660171779821287; // 3/4*sqrt(2); a4
/// let b2  = -0.790569415042094833; // -sqrt(10)/4; -a5
/// let b3  =  0.790569415042094833; // sqrt(10)/4; a5
/// let b4  = -1.060660171779821287; // -3/4*sqrt(2); -a4
/// let b5  =  0.866025403784438647; // sqrt(3)/2; a6
/// let b6  = -0.866025403784438647; // -sqrt(3)/2; -a6
/// let b7  = -0.273861278752583057; // -sqrt(6/5)/4;  -0.25*a7
/// let b8  = -0.612372435695794525; // -sqrt(6)/4;  -0.25*s6
/// let b9  =  1.095445115010332227; // sqrt(6/5); a7
/// let b10 = -0.670820393249936909; // -3/sqrt(5)/2
/// let rel_max = [
/// //  0       1       2     3     4     5     6     7    8     9
/// //  xxx,    xxy,  xxz,  xyy,  xyz,  xzz,  yyy,  yyz, yzz,  zzz  
///     0.0,     b1,  0.0,  0.0,  0.0,  0.0,   b2,  0.0, 0.0,  0.0,  //(3,-3)
///     0.0,    0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0, 0.0,  0.0,  //(3,-2)
///     0.0,     b7,  0.0,  0.0,  0.0,  0.0,   b8,  0.0,  b9,  0.0,  //(3,-1)
///     0.0,    0.0,  b10,  0.0,  0.0,  0.0,  0.0,  b10, 0.0,  1.0,  //(3, 0)
///      b8,    0.0,  0.0,   b7,  0.0,   b9,  0.0,  0.0, 0.0,  0.0,  //(3, 1)
///     0.0,    0.0,   b5,  0.0,  0.0,  0.0,  0.0,   b6, 0.0,  0.0,  //(3, 2)
///      b3,    0.0,  0.0,   b4,  0.0,  0.0,  0.0,  0.0, 0.0,  0.0,  //(3, 3)
/// ].to_vec();
pub const C2S_L3: C2S = C2S::L3(DMatrix10x7 {
    size: [10,7],
    indicing: [1,10],
    data: [
                //  0       1       2     3     4     5     6     7    8     9
        //  xxx,    xxy,  xxz,  xyy,  xyz,  xzz,  yyy,  yyz, yzz,  zzz  
            0.0,     1.060660171779821287,  0.0,  0.0,  0.0,  0.0,   -0.790569415042094833,  0.0, 0.0,  0.0,  //(3,-3)
            0.0,    0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0, 0.0,  0.0,  //(3,-2)
            0.0,     -0.273861278752583057,  0.0,  0.0,  0.0,  0.0,   -0.612372435695794525,  0.0,  1.095445115010332227,  0.0,  //(3,-1)
            0.0,    0.0,  -0.670820393249936909,  0.0,  0.0,  0.0,  0.0,  -0.670820393249936909, 0.0,  1.0,  //(3, 0)
           -0.612372435695794525,    0.0,  0.0,   -0.273861278752583057,  0.0,   1.095445115010332227,  0.0,  0.0, 0.0,  0.0,  //(3, 1)
            0.0,    0.0,   0.866025403784438647,  0.0,  0.0,  0.0,  0.0,   -0.866025403784438647, 0.0,  0.0,  //(3, 2)
            0.790569415042094833,    0.0,  0.0,   -1.060660171779821287,  0.0,  0.0,  0.0,  0.0, 0.0,  0.0,  //(3, 3)
    ]
});

//pub struct DMatrix15x9 {
//    size: [usize;2],
//    indicing: [usize;2],
//    data: [f64;135]
//}

/// The transofrmation matrix from Cartesian to spheric for L=4 (G)
///        let c1  = -0.878310065653679861;  // 3\sqrt(3)/\sqrt(35)
///        let c2  =  0.219577516413419965;  // -a1/4
///        let c3  =  1.195228609334393640;  // sqrt(10/7)
///        let c4  = -0.896421457000795230;  // -3/4*sqrt(10/7)
///        let c5  = -0.400891862868636577;  // -3/4*sqrt(2/7)
///        let c6  =  0.981980506061965716;  //3/2*sqrt(3/7)
///        let c7  = -0.981980506061965716;  //-3/2*sqrt(3/7)
///        let c8  = -0.559016994374947424;  //-sqrt(5)/4
///        let c9  =  0.559016994374947424;  // sqrt(5)/4
///        let c10 =  1.133893419027681682;  // 3/sqrt(7)
///        let c11 = -0.422577127364258289; // -sqrt(5/7)/2
///        let c12 =  0.790569415042094833; // sqrt(10)/4
///        let c13 = -1.060660171779821287; // -3sqrt(2)/4
///        let c14 = -0.790569415042094833; // -sqrt(10)/4
///        let c15 =  1.060660171779821287; // 3sqrt(2)/4
///        let c16 =  0.739509972887452005; // sqrt(35)/8
///        let c17 = -1.299038105676657970; // -3sqrt(3)/4
///        let c18 =  1.118033988749894848; // sqrt(5)/2
///        let c19 = -1.118033988749894848; // -sqrt(5)/2
///        let rel_max = [
///        //  0       1       2     3       4     5       6      7     8       9     10    11    12    13    14
///        //  xxxx,   xxxy,   xxxz, xxyy,   xxyz, xxzz,   xyyy,  xyyz, xyzz,   xzzz  yyyy  yyyz  yyzz  yzzz  zzzz
///            0.0,    c18,    0.0,  0.0,    0.0,  0.0,    c19,   0.0,  0.0,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4,-4)
///            0.0,    0.0,    0.0,  0.0,    c15,  0.0,    0.0,   0.0,  0.0,    0.0,  0.0,  c14,  0.0,  0.0,  0.0, //(4,-3)
///            0.0,    c11,    0.0,  0.0,    0.0,  0.0,    c11,   0.0,  c10,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4,-2)
///            0.0,    0.0,    0.0,  0.0,     c5,  0.0,    0.0,   0.0,  0.0,    0.0,  0.0,   c4,  0.0,   c3,  0.0, //(4,-1)
///          0.375,    0.0,    0.0,   c2,    0.0,   c1,    0.0,   0.0,  0.0,    0.0,0.375,  0.0,   c1,  0.0,  1.0, //(4, 0)
///            0.0,    0.0,     c4,  0.0,    0.0,  0.0,    0.0,    c5,  0.0,     c3,  0.0,  0.0,  0.0,  0.0,  0.0, //(4, 1)
///             c8,    0.0,    0.0,  0.0,    0.0,   c6,    0.0,   0.0,  0.0,    0.0,   c9,  0.0,   c7,  0.0,  0.0, //(4, 2)
///            0.0,    0.0,    c12,  0.0,    0.0,  0.0,    0.0,   c13,  0.0,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4, 3)
///            c16,    0.0,    0.0,  c17,    0.0,  0.0,    0.0,   0.0,  0.0,    0.0,  c16,  0.0,  0.0,  0.0,  0.0, //(4, 4)
///        ].to_vec();
/// 
pub const C2S_L4: C2S = C2S::L4(DMatrix15x9 {
    size: [15,9],
    indicing: [1,15],
    data: [
        //  0       1       2     3       4     5       6      7     8       9     10    11    12    13    14
        //  xxxx,   xxxy,   xxxz, xxyy,   xxyz, xxzz,   xyyy,  xyyz, xyzz,   xzzz  yyyy  yyyz  yyzz  yzzz  zzzz
            0.0,   1.118033988749894848,    0.0,  0.0,    0.0,  0.0,   -1.118033988749894848,   0.0,  0.0,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4,-4)
            0.0,    0.0,    0.0,  0.0,   1.060660171779821287,  0.0,    0.0,   0.0,  0.0,    0.0,  0.0, -0.790569415042094833,  0.0,  0.0,  0.0, //(4,-3)
            0.0,   -0.422577127364258289,    0.0,  0.0,    0.0,  0.0,   -0.422577127364258289,   0.0, 1.133893419027681682,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4,-2)
            0.0,    0.0,    0.0,  0.0,    -0.400891862868636577,  0.0,    0.0,   0.0,  0.0,    0.0,  0.0,  -0.896421457000795230,  0.0,  1.195228609334393640,  0.0, //(4,-1)
          0.375,    0.0,    0.0,   0.219577516413419965,    0.0,   -0.878310065653679861,    0.0,   0.0,  0.0,    0.0,0.375,  0.0,   -0.878310065653679861,  0.0,  1.0, //(4, 0)
            0.0,    0.0,    -0.896421457000795230,  0.0,    0.0,  0.0,    0.0,   -0.400891862868636577,  0.0,    1.195228609334393640,  0.0,  0.0,  0.0,  0.0,  0.0, //(4, 1)
            -0.559016994374947424,    0.0,    0.0,  0.0,    0.0,  0.981980506061965716,    0.0,   0.0,  0.0,    0.0,  0.559016994374947424,  0.0,  -0.981980506061965716,  0.0,  0.0, //(4, 2)
            0.0,    0.0,   0.790569415042094833,  0.0,    0.0,  0.0,    0.0,  -1.060660171779821287,  0.0,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4, 3)
           0.739509972887452005,    0.0,    0.0, -1.299038105676657970,    0.0,  0.0,    0.0,   0.0,  0.0,    0.0, 0.739509972887452005,  0.0,  0.0,  0.0,  0.0, //(4, 4)
    ]
});
/// The transformation matrix from Cartesian to spheric for L=5 (H)
///  let c01 = -1.0910894511799618;  // -5/sqrt(21)
///  let c02 =  0.625;               // 5/8
///  let c03 =  0.36596252735569995; // sqrt(15/7)/4
///  let c11 =  1.2909944487358056;  // sqrt(5/3)
///  let c12 = -1.267731382092775;   //-sqrt(5/7)*3/2
///  let c13 = -0.5669467095138407;  //-sqrt(1/7)*3/2
///  let c14 =  0.4841229182759271;  // sqrt(15)/8
///  let c15 =  0.1613743060919757;  // sqrt(5/3)/8
///  let c16 =  0.21128856368212914; // sqrt(5/7)/4
///  let c21 =  1.118033988749895;   // sqrt(5/4)
///  let c22 =  1.2909944487358056;  // sqrt(5/3)
///  let c23 = -0.8539125638299665;  //-sqrt(35/3)/4
///  let c24 = -0.6454972243679028;  //-sqrt(5/3)/2
///  let c31 =  0.9128709291752769;  // sqrt(10/3)/2
///  let c32 = -1.224744871391589;   //-sqrt(6)/2
///  let c33 = -0.5229125165837972;  //-sqrt(70)/16
///  let c34 =  0.22821773229381923; // sqrt(10/3)/8
///  let c41 =  0.739509972887452;   // sqrt(35)/8
///  let c42 = -1.299038105676658;   //-sqrt(3)*3/4
///  let c43 =  1.118033988749895;   // sqrt(5)/2
///  let c51 =  0.701560760020114;   // sqrt(14)*3/16
///  let c52 =  1.1692679333668567;  // sqrt(14)*5/16
///  let c53 = -1.5309310892394863;  //-sqrt(6)*5/8
///        //   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19   20
///        // 500  410  401  320  311  302  230  221  212  203  140  131  122  113  104  050  041  032  023  014  005
///           0.0, c52, 0.0, 0.0, 0.0, 0.0, c53, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, c51, 0.0, 0.0, 0.0, 0.0, 0.0, //(5,-5)
///           0.0, 0.0, 0.0, 0.0, c43, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-c43, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5,-4)
///           0.0, c33, 0.0, 0.0, 0.0, 0.0,-c34, 0.0,-c32, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-c33, 0.0,-c31, 0.0, 0.0, 0.0, //(5,-3)
///           0.0, 0.0, 0.0, 0.0, c24, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, c24, 0.0, c22, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5,-2)
///           0.0, c15, 0.0, 0.0, 0.0, 0.0, c16, 0.0, c13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, c14, 0.0, c12, 0.0, c11, 0.0, //(5,-1)
///           0.0, 0.0, c02, 0.0, 0.0, 0.0, 0.0, c03, 0.0, c01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, c02, 0.0, c01, 0.0, 1.0, //(5, 0)
///           c14, 0.0, 0.0, c16, 0.0, c12, 0.0, 0.0, 0.0, 0.0, c15, 0.0, c13, 0.0, c11, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5, 1)
///           0.0, 0.0, c23, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, c21, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-c23, 0.0,-c21, 0.0, 0.0, //(5, 2)
///           c33, 0.0, 0.0, c34, 0.0, c31, 0.0, 0.0, 0.0, 0.0,-c33, 0.0, c32, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5, 3)
///           0.0, 0.0, c41, 0.0, 0.0, 0.0, 0.0, c42, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, c41, 0.0, 0.0, 0.0, 0.0, //(5, 4)
///           c51, 0.0, 0.0, c53, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, c52, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5, 5)
/// 
pub const C2S_L5: C2S = C2S::L5(DMatrix21x11 {
    size: [21,11],
    indicing: [1,21],
    data: [
        //   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19   20
        // 500  410  401  320  311  302  230  221  212  203  140  131  122  113  104  050  041  032  023  014  005
           0.0, 1.1692679333668567, 0.0, 0.0, 0.0, 0.0, -1.5309310892394863, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.701560760020114, 0.0, 0.0, 0.0, 0.0, 0.0, //(5,-5)
           0.0, 0.0, 0.0, 0.0, 1.118033988749895, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-1.118033988749895, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5,-4)
           //0.0, -0.5229125165837972, 0.0, 0.0, 0.0, 0.0,-0.22821773229381923, 0.0,0.8660254037844386, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.5229125165837972, 0.0,-0.9128709291752769, 0.0, 0.0, 0.0, //(5,-3)
           0.0, -0.5229125165837972, 0.0, 0.0, 0.0, 0.0,-0.22821773229381923, 0.0,1.224744871391589, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.5229125165837972, 0.0,-0.9128709291752769, 0.0, 0.0, 0.0, //(5,-3)
           0.0, 0.0, 0.0, 0.0, -0.6454972243679028, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.6454972243679028, 0.0, 1.2909944487358056, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5,-2)
           0.0, 0.1613743060919757, 0.0, 0.0, 0.0, 0.0, 0.21128856368212914, 0.0, -0.5669467095138407, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4841229182759271, 0.0, -1.267731382092775, 0.0, 1.2909944487358056, 0.0, //(5,-1)
           0.0, 0.0, 0.625, 0.0, 0.0, 0.0, 0.0, 0.36596252735569995, 0.0, -1.0910894511799618, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.625, 0.0, -1.0910894511799618, 0.0, 1.0, //(5, 0)
           0.4841229182759271, 0.0, 0.0, 0.21128856368212914, 0.0, -1.267731382092775, 0.0, 0.0, 0.0, 0.0, 0.1613743060919757, 0.0, -0.5669467095138407, 0.0, 1.2909944487358056, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5, 1)
           0.0, 0.0, -0.8539125638299665, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.118033988749895, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.8539125638299665, 0.0,-1.118033988749895, 0.0, 0.0, //(5, 2)
           //-0.5229125165837972, 0.0, 0.0, 0.22821773229381923, 0.0, 0.9128709291752769, 0.0, 0.0, 0.0, 0.0,0.5229125165837972, 0.0, -0.8660254037844386, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5, 3)
           -0.5229125165837972, 0.0, 0.0, 0.22821773229381923, 0.0, 0.9128709291752769, 0.0, 0.0, 0.0, 0.0,0.5229125165837972, 0.0, -1.224744871391589, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5, 3)
           0.0, 0.0, 0.739509972887452, 0.0, 0.0, 0.0, 0.0, -1.299038105676658, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.739509972887452, 0.0, 0.0, 0.0, 0.0, //(5, 4)
           0.701560760020114, 0.0, 0.0, -1.5309310892394863, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1692679333668567, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //(5, 5)
    ]
});

impl C2S {
    pub fn to_matrixfullslice(&self) -> MatrixFullSlice<f64> {
        match &self {
            C2S::L0(matr) => {
                MatrixFullSlice {
                    size: &matr.size,
                    indicing: &matr.indicing,
                    data: &matr.data
                }
            },
            C2S::L1(matr) => {
                MatrixFullSlice {
                    size: &matr.size,
                    indicing: &matr.indicing,
                    data: &matr.data
                }
            },
            C2S::L2(matr) => {
                MatrixFullSlice {
                    size: &matr.size,
                    indicing: &matr.indicing,
                    data: &matr.data
                }
            },
            C2S::L3(matr) => {
                MatrixFullSlice {
                    size: &matr.size,
                    indicing: &matr.indicing,
                    data: &matr.data
                }
            },
            C2S::L4(matr) => {
                MatrixFullSlice {
                    size: &matr.size,
                    indicing: &matr.indicing,
                    data: &matr.data
                }
            },
            C2S::L5(matr) => {
                MatrixFullSlice {
                    size: &matr.size,
                    indicing: &matr.indicing,
                    data: &matr.data
                }
            },
        }
    }
}

pub fn c2s_matrix_const(l:usize) -> crate::constants::c2s::C2S {
    if l == 0 {
        C2S_L0
    } else if l ==1 {
        C2S_L1
    } else if l== 2 {
        C2S_L2
    } else if l== 3 {
        C2S_L3
    } else if l== 4 {
        C2S_L4
    } else if l== 5 {
        C2S_L5
    } else {
        panic!("No C2S transformation implementation for l > 5")
    }
}

pub fn c2s_matrix(l:usize) -> MatrixFull<f64> {
    let len_c = (l+1)*(l+2)/2;
    let len_s = 2*l+1;
    if l==0 {
        return MatrixFull::new([len_c,len_s],1.0)
    } else if l==1 {
        let rel_max = [
        //  x    y    z
            1.0, 0.0, 0.0,
            0.0, 1.0, 0.0,
            0.0, 0.0, 1.0
        ].to_vec();
        return unsafe{MatrixFull::from_vec_unchecked([len_c,len_s], rel_max)}
    } else if l==2 {
        let a1 =  0.866025403784438647; //sqrt(3)/2
        let a2 = -0.866025403784438647; //-sqrt(3)/2
        let rel_max = [
        //  0    1    2    3    4    5
        //  x^2, xy,  xz,  y^2, yz,  z^2  
            0.0, 1.0, 0.0, 0.0, 0.0, 0.0,  //(2,-2)
            0.0, 0.0, 0.0, 0.0, 1.0, 0.0,  //(2,-1)
           -0.5, 0.0, 0.0,-0.5, 0.0, 1.0,  //(2, 0)
            0.0, 0.0, 1.0, 0.0, 0.0, 0.0,  //(2, 1)
             a1, 0.0, 0.0,  a2, 0.0, 0.0,  //(2, 2)
            ].to_vec();
        return unsafe{MatrixFull::from_vec_unchecked([len_c,len_s], rel_max)}
    } else if l==3 {
        //let s2 = (2.0f64).sqrt();
        //let s3 = (3.0f64).sqrt();
        //let s5 = (5.0f64).sqrt();
        //let s6 = (6.0f64).sqrt();
        //let a1 = 0.5*s5/s2;  // 1/2\sqrt(5/2)
        //let a2 = 0.5*s3/s2;  // 1/2\sqrt(3/2)
        //let a3 = 0.5*s3*s5;  // 1/2\sqrt(15)
        
        //let a4 = 0.75*s2;
        //let a5 = 0.25*s2*s5;
        //let a6 = 0.5*s3;
        //let a7 = s6/s5;

        //let rel_max = [
        // //  0       1          2         3     4        5         6       7       8       9
        // //  xxx,    xxy,       xxz,      xyy,  xyz,     xzz,      yyy,    yyz,    yzz,    zzz  
        //     0.0,       a4,     0.0,      0.0,  0.0,     0.0,      -a5,    0.0,    0.0,    0.0,  //(3,-3)
        //     0.0,      0.0,     0.0,      0.0,  1.0,     0.0,      0.0,    0.0,    0.0,    0.0,  //(3,-2)
        //     0.0, -0.25*a7,     0.0,      0.0,  0.0,     0.0, -0.25*s6,    0.0,     a7,    0.0,  //(3,-1)
        //     0.0,      0.0, -1.5/s5,      0.0,  0.0,     0.0,      0.0,-1.5/s5,    0.0,    1.0,  //(3, 0)
        //-0.25*s6,      0.0,     0.0, -0.25*a7,  0.0,      a7,      0.0,    0.0,    0.0,    0.0,  //(3, 1)
        //     0.0,      0.0,      a6,      0.0,  0.0,     0.0,      0.0,    -a6,    0.0,    0.0,  //(3, 2)
        //      a5,      0.0,     0.0,      -a4,  0.0,     0.0,      0.0,    0.0,    0.0,    0.0,  //(3, 3)
        //].to_vec();

        let b1  =  1.060660171779821287; // 3/4*sqrt(2); a4
        let b2  = -0.790569415042094833; // -sqrt(10)/4; -a5
        let b3  =  0.790569415042094833; // sqrt(10)/4; a5
        let b4  = -1.060660171779821287; // -3/4*sqrt(2); -a4
        let b5  =  0.866025403784438647; // sqrt(3)/2; a6
        let b6  = -0.866025403784438647; // -sqrt(3)/2; -a6
        let b7  = -0.273861278752583057; // -sqrt(6/5)/4;  -0.25*a7
        let b8  = -0.612372435695794525; // -sqrt(6)/4;  -0.25*s6
        let b9  =  1.095445115010332227; // sqrt(6/5); a7
        let b10 = -0.670820393249936909; // -3/sqrt(5)/2
        let rel_max = [
        //  0       1       2     3     4     5     6     7    8     9
        //  xxx,    xxy,  xxz,  xyy,  xyz,  xzz,  yyy,  yyz, yzz,  zzz  
            0.0,     b1,  0.0,  0.0,  0.0,  0.0,   b2,  0.0, 0.0,  0.0,  //(3,-3)
            0.0,    0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0, 0.0,  0.0,  //(3,-2)
            0.0,     b7,  0.0,  0.0,  0.0,  0.0,   b8,  0.0,  b9,  0.0,  //(3,-1)
            0.0,    0.0,  b10,  0.0,  0.0,  0.0,  0.0,  b10, 0.0,  1.0,  //(3, 0)
             b8,    0.0,  0.0,   b7,  0.0,   b9,  0.0,  0.0, 0.0,  0.0,  //(3, 1)
            0.0,    0.0,   b5,  0.0,  0.0,  0.0,  0.0,   b6, 0.0,  0.0,  //(3, 2)
             b3,    0.0,  0.0,   b4,  0.0,  0.0,  0.0,  0.0, 0.0,  0.0,  //(3, 3)
        ].to_vec();
        return unsafe{MatrixFull::from_vec_unchecked([len_c,len_s], rel_max)}
    } else if l==4 {
        let a1  = -0.878310065653679861;  // 3\sqrt(3)/\sqrt(35)
        let a2  =  0.219577516413419965;  // -a1/4
        let a3  =  1.195228609334393640;  // sqrt(10/7)
        let a4  = -0.896421457000795230;  // -3/4*sqrt(10/7)
        let a5  = -0.400891862868636577;  // -3/4*sqrt(2/7)
        let a6  =  0.981980506061965716;  //3/2*sqrt(3/7)
        let a7  = -0.981980506061965716;  //-3/2*sqrt(3/7)
        let a8  = -0.559016994374947424;  //-sqrt(5)/4
        let a9  =  0.559016994374947424;  // sqrt(5)/4
        let a10 =  1.133893419027681682;  // 3/sqrt(7)
        let a11 = -0.422577127364258289; // -sqrt(5/7)/2
        let a12 =  0.790569415042094833; // sqrt(10)/4
        let a13 = -1.060660171779821287; // -3sqrt(2)/4
        let a14 = -0.790569415042094833; // -sqrt(10)/4
        let a15 =  1.060660171779821287; // 3sqrt(2)/4
        let a16 =  0.739509972887452005; // sqrt(35)/8
        let a17 = -1.299038105676657970; // -3sqrt(3)/4
        let a18 =  1.118033988749894848; // sqrt(5)/2
        let a19 = -1.118033988749894848; // -sqrt(5)/2
        let rel_max = [
        //  0       1       2     3       4     5       6      7     8       9     10    11    12    13    14
        //  xxxx,   xxxy,   xxxz, xxyy,   xxyz, xxzz,   xyyy,  xyyz, xyzz,   xzzz  yyyy  yyyz  yyzz  yzzz  zzzz
            0.0,    a18,    0.0,  0.0,    0.0,  0.0,    a19,   0.0,  0.0,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4,-4)
            0.0,    0.0,    0.0,  0.0,    a15,  0.0,    0.0,   0.0,  0.0,    0.0,  0.0,  a14,  0.0,  0.0,  0.0, //(4,-3)
            0.0,    a11,    0.0,  0.0,    0.0,  0.0,    a11,   0.0,  a10,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4,-2)
            0.0,    0.0,    0.0,  0.0,     a5,  0.0,    0.0,   0.0,  0.0,    0.0,  0.0,   a4,  0.0,   a3,  0.0, //(4,-1)
          0.375,    0.0,    0.0,   a2,    0.0,   a1,    0.0,   0.0,  0.0,    0.0,0.375,  0.0,   a1,  0.0,  1.0, //(4, 0)
            0.0,    0.0,     a4,  0.0,    0.0,  0.0,    0.0,    a5,  0.0,     a3,  0.0,  0.0,  0.0,  0.0,  0.0, //(4, 1)
             a8,    0.0,    0.0,  0.0,    0.0,   a6,    0.0,   0.0,  0.0,    0.0,   a9,  0.0,   a7,  0.0,  0.0, //(4, 2)
            0.0,    0.0,    a12,  0.0,    0.0,  0.0,    0.0,   a13,  0.0,    0.0,  0.0,  0.0,  0.0,  0.0,  0.0, //(4, 3)
            a16,    0.0,    0.0,  a17,    0.0,  0.0,    0.0,   0.0,  0.0,    0.0,  a16,  0.0,  0.0,  0.0,  0.0, //(4, 4)
        ].to_vec();
        return unsafe{MatrixFull::from_vec_unchecked([len_c,len_s], rel_max)}

    } else {
        panic!("the angular momentum larger than 4 is not yet implemented for DFA");
        return MatrixFull::empty()
    }

}